#Base plate pull-up model based on Malhotra and Veletsos (1994) as implemented in Bakalis et al. (2017)
 
import openseespy.opensees as op
import vfo.vfo as vfo
import os
import numpy as np
import matplotlib.pyplot as plt
import time
import Units as ut

op.wipe();

print('Model generation started...')

direc = './ResultsBP'

if(os.path.exists(direc)=='False'):
    os.mkdir(direc);


R = 13.9*ut.m #Tank Radius in m
H = 16.5*ut.m #Contained Liquid Height in m
n = 8 #Number of beams used to modelling the base plate, Bakalis et al. (2017)
tw = 17.7*ut.mm #wall thickness in mm
tb = 6.4*ut.mm #base plate thickness
ta = 9.5*ut.mm #annular ring thickness
ro = 998*ut.kg_m3 #liquid density in kg/m3  
Ew = 1000*ut.MPa #Young Modulus of rigid foundation in MPa 
Es = 210000*ut.MPa #Young Modulus of steel in MPa
Fys = 235*ut.MPa #Yield strength of steel, MPa
epy = Fys/Es
Bs = 0.01 #strain hardening ratio
v = 0.3 #Poisson ratio of steel

fillr = 0.95 #Pencentage of tank's height filled with liquid 

bw = 2*np.pi*R/n  #uniform width, Bakalis et al. (2017)
ros = 7850*ut.kg_m3 #density of steel
Ww = ros*ut.g*H*bw*tw #weight of wall over one side of the strip
mr = 35*ut.ton
Wr = mr*ut.g #weight of roof
Wrr = Wr/n
ph = ro*ut.g*H*fillr  #Hydrostatic pressure at tank bottom

ktt = Es*bw*tw**2*(tw/R)**0.5/(2*(3*(1-v**2))**0.75) #Rotational stiffness due to wall, Bakalis et al. (2017)
kuu = Es*bw*(tw/R)**1.5/((3*(1-v**2))**0.25)  #Axial stiffness due to wall, Bakalis et al. (2017) 
Mr = bw*R*tw*ph/(2*(3*(1-v**2))**0.5)   #Point moment at bas plate wall junction due to hydrostatic pressure, Bakalis et al. (2017)
Nr = bw*(R*tw)**0.5*ph/(2*(3*(1-v**2))**0.25) #Point radial load due to hydrostatic pressure, Bakalis et al. (2017) 

kw = Ew*bw #stiffness of zero-length Winkler spring at tank foundation interface

lb = 50*ut.mm  #length of beam sub-element aprox. 15tb, Bakalis et al. (2017) 
nb = int(2*R/lb)   #number of beam sub-elements

la = 2000*ut.mm    #length of annular ring 
na = int(round(la/lb,0))   #number of elements in the annular rings 
wph = ph*bw  #distributed load on beam due to hydrostatic pressure
print(wph)
#Total Weight of Liquid
Wtot = wph*2*R
# Point Load in Nodes
wpn = Wtot/(nb+1)


op.model('basic', '-ndm', 2, '-ndf', 3) 

#Define nodes
nodeTag1 = 10000
nodeTag2 = 20000

for i in range(nb+1):
    op.node(nodeTag1+i,i*lb,0) 
 

#op.fix(nodeTag1+nb,1,0,1)    

for i in range(1,nb):
    op.node(nodeTag2+i,i*lb,0)
    op.fix(nodeTag2+i,1,1,1)
    
    
op.fix(nodeTag1,1,1,0) #pined support at compression edge

nodeMaster = 30000

op.node(nodeMaster,nb*lb,0)
op.fix(nodeMaster,1,0,1)

#Define Materials
steelMatTag = 1
op.uniaxialMaterial('Steel01', steelMatTag, Fys, Es, Bs)
winkMatTag = 2
op.uniaxialMaterial('ENT',winkMatTag,kw)
kttMatTag = 3
op.uniaxialMaterial('Elastic',kttMatTag,ktt)
kuuMatTag = 4
op.uniaxialMaterial('Elastic',kuuMatTag,kuu)
rigidMatTag = 5
op.uniaxialMaterial('Elastic',rigidMatTag,ut.Ubig)
flexMatTag = 6
op.uniaxialMaterial('Elastic',flexMatTag,ut.Usmall)

#Create base plate fiber section
secTag = 1
op.section('Fiber',secTag)
# number of fibers in z-direction
nfZ = 4
# number of fibers in y-direction
nfY = int(round(bw/5,0))
y1 = -bw/2
yJ = bw/2
z1 = -tb/2
zJ = tb/2
op.patch('rect', steelMatTag, nfZ, nfY, z1, y1, zJ,yJ)
#Create annular ring fiber section
secTaga = 2
op.section('Fiber',secTaga)
# number of fibers in z-direction
y1a = -bw/2
yJa = bw/2
z1a = -ta/2
zJa = ta/2
op.patch('rect', steelMatTag, nfZ, nfY, z1a, y1a, zJa,yJa)

#Create elements

#Winkler Springs
eleTagW = 20000
for i in range(1,nb):
    op.element('zeroLength',eleTagW+(i-1),nodeTag2+i,nodeTag1+i,'-mat',winkMatTag,'-dir',2)

#Beam Elements
intgTag = 1
intgTaga = 2
numInt = 10
maxIter = 10
tol = 1e-12
#Define geometric transformation
BPTransf = 'Corotational'
transfTag = 1
op.geomTransf(BPTransf,transfTag)

op.beamIntegration('Legendre',intgTag,secTag,numInt)
op.beamIntegration('Legendre',intgTaga,secTaga,numInt)
eleTagBeam = 30000
for i in range(0,nb):
    if (i<na or i>=nb-na):
        op.element('forceBeamColumn',eleTagBeam+i,nodeTag1+i,nodeTag1+(i+1),transfTag,intgTaga,'-iter',maxIter,tol)
    else:
        op.element('forceBeamColumn',eleTagBeam+i,nodeTag1+i,nodeTag1+(i+1),transfTag,intgTag,'-iter',maxIter,tol)
ktueleTag = 1
op.element('zeroLength',ktueleTag,nodeMaster,nodeTag1+nb,'-mat',kuuMatTag,rigidMatTag,kttMatTag,'-dir',1,2,3)


#Run Gravity Analysis
constTS = 1
op.timeSeries('Constant', constTS)
patternTagG = 1
op.pattern('Plain',patternTagG,constTS)

for i in range(nb+1):
      op.load(nodeTag1+i,0,-wpn,0)
#op.eleLoad('-range',eleTagBeam,eleTagBeam+nb-1,'-type','-beamUniform',-wph)
op.load(nodeTag1+nb,Nr,0,-Mr)
#op.load(nodeTag1+nb,0,-(Ww+Wrr),0)

Tol = 1e-8 # convergence tolerance for test
NstepGravity = 10
DGravity = 1/NstepGravity
op.integrator('LoadControl', DGravity) # determine the next time step for an analysis
op.numberer('RCM') # renumber dof's to minimize band-width (optimization), if you want to
op.system('SparseGeneral', '-piv') # how to store and solve the system of equations in the analysis
op.constraints('Plain') # how it handles boundary conditions
op.test('EnergyIncr', Tol, 50) # determine if convergence has been achieved at the end of an iteration step
op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
op.analysis('Static') # define type of analysis static or transient
op.analyze(NstepGravity) # apply gravity
op.loadConst('-time', 0.0) #maintain constant gravity loads and reset time to zero
print('Static analysis completed')


#Define recorders
op.recorder('Node', '-file', direc+'/disp.out','-time', '-node', nodeMaster, '-dof', 2, 'disp')
op.recorder('Node', '-file', direc+'/WinklerSprings.out','-time', '-nodeRange', nodeTag2+1, nodeTag2+nb-1,'-dof', 2, 'reaction')
op.recorder('Node', '-file', direc+'/WinklerSprings_disp.out','-time', '-nodeRange', nodeTag2+1, nodeTag2+nb-1,'-dof', 2, 'disp')


#Run Pullup Analysis

IDctrlNode = nodeMaster			# node where disp is read for disp control
IDctrlDOF  = 2			        # degree of freedom read for disp control (1 = x displacement)
Dmax  = 1000*ut.mm	            # maximum displacement of pushover
Dincr = Dmax/10000  # displacement increment
Tol = 1e-8 # convergence tolerance for test
print('Running Pushover...')


POtag = 100;
linTS = 2
op.timeSeries('Linear',linTS)

op.pattern('Plain', POtag, linTS)     
op.load(IDctrlNode,0,1,0)


op.constraints('Plain')
op.numberer('RCM')
op.system('BandGeneral')
op.test('EnergyIncr', Tol, 50)
op.algorithm('Newton')
op.integrator('DisplacementControl', IDctrlNode, IDctrlDOF, Dincr)
op.analysis('Static')

ok = 0
currentDisp = 0.0
Nsteps = int(Dmax/Dincr)		# number of pushover analysis steps
ok = op.analyze(Nsteps)        # this will return zero if no convergence problems were encountered
# ----------------------------------------------if convergence failure-------------------------
if(ok != 0):   
    # if analysis fails, we try some other stuff, performance is slower inside this loop
    Dstep =  0.0
    ok = 0
    while(Dstep <= 1.0 and ok == 0):
        controlDisp = op.nodeDisp(IDctrlNode,IDctrlDOF)
        Dstep = controlDisp/Dmax
        ok = op.analyze(1)
        # if analysis fails, we try some other stuff
        # performance is slower inside this loop	global maxNumIterStatic;	    # max no. of iterations performed before "failure to converge" is ret'd
        if(ok != 0): 
            print("Trying Newton with Initial Tangent ..")
            op.test('NormDispIncr',1.0e-5, 20000, 0)
            op.algorithm('Newton', '-initial')
            ok = op.analyze(1)
            op.test('EnergyIncr', Tol, 100)
            op.algorithm('Newton')
            print("Back to original parameters ..")           
        if(ok != 0):
            print("Trying Newton with Line Search ..")
            op.test('RelativeEnergyIncr', 10e-6, 1000) 
            op.algorithm('NewtonLineSearch')
            ok = op.analyze(1)
            op.test('EnergyIncr', 1.0e-10, 100)
            op.algorithm('Newton')
            print("Back to original parameters ..")
        if(ok != 0):
            print("Trying Broyden ..")
            op.algorithm('Broyden', 8, 100)
            ok = op.analyze(1)
            op.algorithm('Newton')
            print("Back to original parameters ..")
        if(ok != 0): 
            print("Trying NewtonWithLineSearch ..")
            op.algorithm('NewtonLineSearch', 10e-4, 10000)
            ok = op.analyze(1)
            op.algorithm('Newton')
            print("Back to original parameters ..")
        # end while loop
    # end if ok !0

if(ok != 0):
	print('Problem with Pushover')
else:
	print('Done')
		
#vfo.createODB(model="BasePlate",loadcase="Pushover")
#vfo.plot_deformedshape(model="BasePlate",loadcase="Pushover")
#opv.plot_model(node_labels=0,element_labels=0)
#vfo.plot_model(show_nodes = "yes")

#plt.show()





